#ifndef neuread_H
#define neuread_H

#include <fstream>
#include <map>
#include <vector>

#include "geometry.h"
#include "mesh.h"
#include "boundary.h"

namespace ABAN
{

    class EFileNotFound : public EInputErr{};
    class EForbiddenElementType : public EInputErr{};

    //this class reads gambit neutral file in combination with an ABANBOUNDARY file for boundary conditions
    class ReadNEU
    {
        private:
        struct BCTypeData
        {
            SeqType type;
            std::vector<FloatType> values;
        };
        std::ifstream neuFile,bcFile;
        SeqType nodeNum,cellNum,grpNum,bcNum/*boundary condition*/;
        std::map<std::string,BCTypeData> bcValues;
        //this container is just used for determining faces of boundary conditions
        DatumContainer<CellConnectivity> nonSortedConnectivity;
        //this container brings a direct access to cells for us
        //(cos aCells is list, we can not access its elements by operator[])
        vector<CellIterator> cellsDirectAccess;

        void bypass(SeqType n=1);
        void readBC();

        void readNEUHeader();
        void readNEUNodes(VectorContainer *aNodes,CellContainer *aCells);
        void readNEUCells(VectorContainer *aNodes,CellContainer *aCells);
        void bypassNEUGroupSection();
        void readNEUBoundaries(VectorContainer *aNodes,CellContainer *aCells);

        void readWholeNEU(VectorContainer *aNodes,CellContainer *aCells);

        public:
        ReadNEU(std::string neuFileName,std::string bcFileName);
        ~ReadNEU();

        static void run(VectorContainer *aNodes,CellContainer *aCells,std::string neuFileName,std::string bcFileName);
    };
    //-------------------------------------------------------
    ReadNEU::ReadNEU(std::string neuFileName,std::string bcFileName)
    {
        neuFile.open(neuFileName.c_str());
        bcFile.open(bcFileName.c_str());
    }
    //-------------------------------------------------------
    ReadNEU::~ReadNEU()
    {
        bcFile.close();
        neuFile.close();
    }
    //-------------------------------------------------------
    void ReadNEU::bypass(SeqType n)
    {
        std::string dummyStr;
        for(SeqType i=0;i<n;i++) neuFile>>dummyStr;
    }
    //-------------------------------------------------------
    void ReadNEU::readBC()
    {
        SeqType noOfBCValues;
        bcFile>>noOfBCValues;
        for (SeqType boundaryCounter=0;boundaryCounter<noOfBCValues;boundaryCounter++)
        {
            std::string bcName;
            bcFile>>bcName;
            BCTypeData bcData;
            bcFile>>bcData.type;
            SeqType noOfVecElems;
            bcFile>>noOfVecElems;
            FloatType parameter;
            for (SeqType vecElemCounter=0;vecElemCounter<noOfVecElems;vecElemCounter++)
                {
                    bcFile>>parameter;
                    bcData.values.push_back(parameter);
                }
            bcValues[bcName]=bcData;
        }
    }
    //-------------------------------------------------------
    void ReadNEU::readNEUHeader()
    {
        neuFile>>nodeNum;
        neuFile>>cellNum;
        neuFile>>grpNum;
        neuFile>>bcNum;
    }
    //-------------------------------------------------------
    void ReadNEU::readNEUNodes(VectorContainer *aNodes,CellContainer *aCells)
    {
        Vector vec;
        for (int i=0;i<nodeNum;i++)
        {
            bypass();
            for (SeqType j=dirXP;j<noOfAxis;j++) neuFile>>vec[j];
            aNodes->push_back(vec);
        }
    }
    //-------------------------------------------------------
    void ReadNEU::readNEUCells(VectorContainer *aNodes,CellContainer *aCells)
    {
        for (SeqType i=0;i<cellNum;i++)
        {
            SeqType cellNo;
            CellConnectivity gCell;

            neuFile>>cellNo;//Read cell number
            cellNo--;//In opposite to GAMBIT's cell numbering, we start cell numbers from 0
            SeqType cellType;
            neuFile>>cellType;//Read code of cell type
            if (cellType!=4) throw new EForbiddenElementType();
            SeqType nodesPerThisCell;
            neuFile>>nodesPerThisCell;//Read number of nodes per cell
            if (nodesPerThisCell!=nodesPerCell) throw new EForbiddenElementType();

            //Reading cells connectivity
            VectorIterator vi;
            for (int j=0;j<nodesPerCell;j++)
            {
                SeqType nodeNumber;
                neuFile>>nodeNumber;//Read node number
                vi = aNodes->begin() + nodeNumber - 1;
                gCell[j]=vi;
            }

            //push it back before sorting
            nonSortedConnectivity.push_back(gCell);

            //Re-sort nodes in ABAN order
            gCell.rearrangeNodes();

            //Adding a non-adapted cell
            aCells->appendGrandCell(cellNo,gCell);

            CellIterator aCellsLastIterator=aCells->end();
            aCellsLastIterator--;
            cellsDirectAccess.push_back(aCellsLastIterator);
        }// of reading cell data
    }
    //-------------------------------------------------------
    void ReadNEU::bypassNEUGroupSection()
    {
        std::string dummyStr;
        for (SeqType grpCounter=0;grpCounter<grpNum;grpCounter++)
        {
            dummyStr="";
            while (dummyStr!="ENDOFSECTION")
                neuFile>>dummyStr;
        }
    }
    //-------------------------------------------------------
    void ReadNEU::readNEUBoundaries(VectorContainer *aNodes,CellContainer *aCells)
    {
        for (SeqType bcCounter=0;bcCounter<bcNum;bcCounter++)
        {
            SeqType elemsOfThisBC;
            std::string nameOfBC;
            const SeqType gambitBrickNodes [facesPerCell][nodesPerFace]=
              {
                  {0,1,5,4},
                  {1,3,7,5},
                  {3,2,6,7},
                  {2,0,4,6},
                  {1,0,2,3},
                  {4,5,7,6}
              };

            //section header ra dur mirizimmmmmmmm!
            bypass(3);

            neuFile>>nameOfBC;
            bypass();
            neuFile>>elemsOfThisBC;
            bypass();//tedad ghalam dadehaye ezafe dar har radif
            bypass();//neuFile>>typeBC;

            for (SeqType elem=0;elem<elemsOfThisBC;elem++)
            {
                SeqType bcCellNo,bcFaceNo,dir;

                neuFile>>bcCellNo;bcCellNo--;
                bypass();//for this project this value is always 4 (Brick)
                neuFile>>bcFaceNo;bcFaceNo--;

                //converting gambit's face no to ABAN standard no.

                //rtrievive face center from two nodes(0 & 2) on a diameter of face
                Vector faceCenter;
                faceCenter[dirXP]=((*nonSortedConnectivity[bcCellNo][gambitBrickNodes[bcFaceNo][0]])[dirXP]+
                                   (*nonSortedConnectivity[bcCellNo][gambitBrickNodes[bcFaceNo][2]])[dirXP] )/2.0;
                faceCenter[dirYP]=((*nonSortedConnectivity[bcCellNo][gambitBrickNodes[bcFaceNo][0]])[dirYP]+
                                   (*nonSortedConnectivity[bcCellNo][gambitBrickNodes[bcFaceNo][2]])[dirYP] )/2.0;
                faceCenter[dirZP]=((*nonSortedConnectivity[bcCellNo][gambitBrickNodes[bcFaceNo][0]])[dirZP]+
                                   (*nonSortedConnectivity[bcCellNo][gambitBrickNodes[bcFaceNo][2]])[dirZP] )/2.0;

                Vector cellCenter = cellsDirectAccess[bcCellNo]->center;

                if (fabs(cellCenter[dirXP]-faceCenter[dirXP])>zeroRadious) dir=dirXP; else
                if (fabs(cellCenter[dirYP]-faceCenter[dirYP])>zeroRadious) dir=dirYP; else
                if (fabs(cellCenter[dirZP]-faceCenter[dirZP])>zeroRadious) dir=dirZP;

                if (cellCenter[dir]>faceCenter[dir]) dir=oppositeSide[dir];

                //setting boundary condition class
                AbstractBoundaryCondition* bc;
                switch (bcValues[nameOfBC].type)//(typeBC)
                {
                    case 0: //constant velocities and zero gradient neuman prssure
                    {
                        Vector constVels(bcValues[nameOfBC].values[dirXP],
                                         bcValues[nameOfBC].values[dirYP],
                                         bcValues[nameOfBC].values[dirZP]);

                        bc = new ConstVelBoundary(cellsDirectAccess[bcCellNo],dir,constVels);
                        break;
                    }
                    case 1: //simple outlet
                    {
                        bc = new SimpleOutlet(cellsDirectAccess[bcCellNo],dir);
                        break;
                    }

                    case 2: //simple convectcive outlet
                    {
                        bc = new SimpleConvectiveOutlet(cellsDirectAccess[bcCellNo],dir,bcValues[nameOfBC].values[0]);
                        break;
                    }
                }

                //add the boundary condition to the set of boundary conditions
                aCells->boundary->push_back(bc);
            }
            bypass();//endofsection
        }
    }
    //-------------------------------------------------------
    void ReadNEU::readWholeNEU(VectorContainer *aNodes,CellContainer *aCells)
    {
        //Khenzer penzeraye avvalesh ro mirizim door
        bypass(20);
        //Reading Node Number and Cell Number
        readNEUHeader();
        //Khenzer penzeraye vasatesh ro mirizim door
        bypass(6);
        //Reading Nodes' Coordinations
        readNEUNodes(aNodes,aCells);
        //Khenzer penzeraye vasatesh ro mirizim door
        bypass(3);
        //Reading Cell Data
        readNEUCells(aNodes,aCells);
        bypass(); //end of section
        //deleting cell-group data
        bypassNEUGroupSection();
        //reading boundary conditions
        readNEUBoundaries(aNodes,aCells);
    }
    //-------------------------------------------------------
    void ReadNEU::run(VectorContainer *aNodes,CellContainer *aCells,std::string neuFileName,std::string bcFileName)
    {
        ReadNEU readerClass(neuFileName,bcFileName);
        readerClass.readBC();
        readerClass.readWholeNEU(aNodes,aCells);
    }
    //-------------------------------------------------------
}
#endif
